QC and filter seurat object

library_load <- suppressMessages(
    list(
        # Seurat 
        library(Seurat, quietly=TRUE),

        # Single R
        library(SingleCellExperiment, quietly=TRUE),
        library(SingleR, quietly=TRUE),
        library(celldex, quietly=TRUE),

        # Data 
        library(dplyr, quietly=TRUE),
        library(reshape2, quietly=TRUE),
        library(Matrix, quietly=TRUE),

        # Source 
        library(msigdbr, quietly=TRUE),

        # Plotting
        library(ggplot2, quietly=TRUE),
        library(ComplexHeatmap, quietly=TRUE),
        library(gridExtra, quietly=TRUE),
        library(patchwork, quietly=TRUE),
        library(RColorBrewer, quietly=TRUE),
        library(viridis, quietly=TRUE)
    )
)

# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")

# Source files
source("plotting_global.R")
source("bin/seurat_function.R")

# Cell cycle
source("bin/ccGenes.R")
source("bin/CellCycle.R")

# SingleRPlot
source("bin/SingleRPlot.R")

Warning message:
“For function ‘ToCorrMatrix’, signature ‘DFrame’: argument in method definition changed from (corr) to (object)”

options(warn=-1)
ht_opt$message=FALSE

Parameter settings

# Filtering Parameter
nFeature_RNA_min_m <- 300
nFeature_RNA_min_p <- 300

nFeature_RNA_max_m <- 3000
nFeature_RNA_max_p <- 3000

nCount_RNA_min_m <- 1000
nCount_RNA_min_p <- 1000

nCount_RNA_max_m <- 20000
nCount_RNA_max_p <- 20000

pMt_RNA_max_m <- 5
pMt_RNA_max_p <- 7.5

# Files 
so_raw_file <- "data/object/raw.rds"

so_qc_file <- "data/object/seurat.rds"

# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()

Import Seurat object

so_raw <- readRDS(so_raw_file)

Add filter meta data to Seurat object

so_raw$nFeature_RNA_max <- ifelse(so_raw$tissue == "Myeloid", nFeature_RNA_max_m, nFeature_RNA_max_p)
so_raw$nFeature_RNA_min <- ifelse(so_raw$tissue == "Myeloid", nFeature_RNA_min_m, nFeature_RNA_min_p)

so_raw$nCount_RNA_max <- ifelse(so_raw$tissue == "Myeloid", nCount_RNA_max_m, nCount_RNA_max_p)
so_raw$nCount_RNA_min <- ifelse(so_raw$tissue == "Myeloid", nCount_RNA_min_m, nCount_RNA_min_p)

so_raw$pMt_RNA_max <- ifelse(so_raw$tissue == "Myeloid", pMt_RNA_max_m, pMt_RNA_max_p)

# QC 
so_raw$qc_class <- ifelse(
  so_raw$cellranger_class == "Cell" &
    so_raw$nFeature_RNA <= so_raw$nFeature_RNA_max & 
    so_raw$nFeature_RNA >= so_raw$nFeature_RNA_min & 
    so_raw$nCount_RNA <= so_raw$nCount_RNA_max & 
    so_raw$nCount_RNA >= so_raw$nCount_RNA_min &
    so_raw$pMt_RNA <= so_raw$pMt_RNA_max, 
  "pass", "fail"
  )

so_raw$tissue <- factor(so_raw$tissue, levels=c("Myeloid", "Progenitor"))
so_raw$treatment <- factor(so_raw$treatment, levels=c("NaCl", "CpG"))

Filter cells by rank plot

Empty droplets were determined with CellRanger V3.0.2 Lun et al., 2019 EmptyDrop heuristic. RNAse activity of granulocytes might be wrongly identified as empty cells by CellRanger.

Typical Sample A steep drop-off is indicative of good separation between the cell-associated barcodes and the barcodes associated with empty GEMs. A ideal barcode rank plot has a distincitve shape, which is referred to as a "cliff and knee".
Heterogeneous Sample Heterogeneous populations of cells in a sample result in two "cliff and knee" distributions. However, there should still be clear separation between the bacodes.
Compromised Sample Round curve and lack of steep cliff may indicate low sample quality or loss of single-cell behavior. This can be due to a wetting failure, premature cell lysis, or low cell viability.
Compromised Sample Defined cliff and knee, but the total number of barcodes detected may be lower than expected. This can be caused by a sample clog or inaccurate cell count.

options(repr.plot.width=5, repr.plot.height=5)
rank_plot <- ggplot(so_raw@meta.data, aes(x=log10(nCount_RNA_rank), y=log10(nCount_RNA), color=cellranger_class)) + 
  geom_point() + 
  scale_color_manual(values=color$cellranger_class) +
  ggtitle("Barcode rank plot") +
  xlab("log10(cell barcode rank)") + ylab("log10(cell UMI counts)") + 
  facet_grid(tissue~treatment) + 
  theme(aspect.ratio=1, legend.position="bottom")
rank_plot
ggsave(rank_plot, filename="result/plot/seurat/rank_plot.png", width=5, height=5)

Filter by CellRanger class

so_qc <- subset(so_raw, subset=cellranger_class == "Cell")

Filter by UMI and Feature count

The Progenitor NaCl sample show a high MT%. It would be interessting to see whether those are apoptoitic cells or if they have a increased energy demand. That could be done by correlating the MT% with complementary nuclear mitochondrial gene expression from the OXPHOS system.

qc_1 <- ggplot(so_qc@meta.data, aes(x=log10(nCount_RNA), fill=tissue)) + 
  geom_density() + 
  ggtitle("Density plot UMI count") + xlab("log10(UMI count)") + ylab("Density") +
  geom_vline(aes(xintercept=log10(nCount_RNA_min)), color="red", linetype="longdash") +
  geom_vline(aes(xintercept=log10(nCount_RNA_max)), color="red", linetype="longdash") + 
  scale_x_continuous(breaks=integer_breaks()) + 
  scale_fill_manual(values=color$tissue) +
  facet_grid(tissue~treatment+sample_rep) + 
  theme(legend.position="bottom", aspect.ratio=1) 

qc_2 <- ggplot(so_qc@meta.data, aes(x=log10(nFeature_RNA), fill=tissue)) + 
  geom_density() + 
  ggtitle("Density plot Feature count") + xlab("log10(Feature count)") + ylab("Density") +
  geom_vline(aes(xintercept=log10(nFeature_RNA_min)), color="red", linetype="longdash") +
  geom_vline(aes(xintercept=log10(nFeature_RNA_max)), color="red", linetype="longdash") + 
  scale_x_continuous(breaks=integer_breaks()) + 
  scale_fill_manual(values=color$tissue) +
  facet_grid(tissue~treatment+sample_rep) + 
  theme(legend.position="bottom", aspect.ratio=1)

qc_3 <- ggplot(so_qc@meta.data, aes(x=pMt_RNA, fill=tissue)) + 
  geom_density() + 
  ggtitle("Density plot Mt %") + xlab("Mt [%]") + ylab("Density") +
  geom_vline(aes(xintercept=pMt_RNA_max), color="red", linetype="longdash") +
  scale_x_continuous(breaks=integer_breaks()) +
  scale_fill_manual(values=color$tissue) +
  facet_grid(tissue~treatment+sample_rep)  + 
  theme(legend.position="bottom", aspect.ratio=1)

options(repr.plot.width=15, repr.plot.height=10)
qc_1 + qc_2 + qc_3 + plot_layout(ncol=2) & theme(legend.position="bottom")
ggsave(qc_1, filename="result/plot/seurat/density_umi.png", width=5, height=5)
ggsave(qc_2, filename="result/plot/seurat/density_feature.png", width=5, height=5)
ggsave(qc_3, filename="result/plot/seurat/density_mt.png", width=5, height=5)

sc_1 <- ggplot(so_qc@meta.data, aes(x=log10(nCount_RNA), y=log10(nFeature_RNA), color=pMt_RNA)) + 
  geom_point() + ggtitle("Mitochondrial gene percentage") + ylab("log10(feature count)") + xlab("log10(umi count)") + 
  geom_vline(aes(xintercept=log10(nCount_RNA_min)), color="red", linetype="longdash") +
  geom_vline(aes(xintercept=log10(nCount_RNA_max)), color="red", linetype="longdash") +
  geom_hline(aes(yintercept=log10(nFeature_RNA_min)), color="red", linetype="longdash") + 
  geom_hline(aes(yintercept=log10(nFeature_RNA_max)), color="red", linetype="longdash") + 
  facet_grid(tissue~treatment+sample_rep) + theme(aspect.ratio=1, legend.position="bottom") + 
  scale_size(guide=guide_legend(direction="vertical"))

sc_2 <- ggplot(so_qc@meta.data, aes(x=log10(nCount_RNA), y=log10(nFeature_RNA), color=pHb_RNA)) + 
  geom_point() + ggtitle("Hemoglobin gene percentage") + ylab("log10(feature count)") + xlab("log10(umi count)") +
  geom_vline(aes(xintercept=log10(nCount_RNA_min)), color="red", linetype="longdash") +
  geom_vline(aes(xintercept=log10(nCount_RNA_max)), color="red", linetype="longdash") +
  geom_hline(aes(yintercept=log10(nFeature_RNA_min)), color="red", linetype="longdash") + 
  geom_hline(aes(yintercept=log10(nFeature_RNA_max)), color="red", linetype="longdash") + 
  facet_grid(tissue~treatment+sample_rep) + theme(aspect.ratio=1, legend.position="bottom") + 
  scale_size(guide=guide_legend(direction="vertical"))

sc_3 <- ggplot(so_qc@meta.data, aes(x=log10(nCount_RNA), y=log10(nFeature_RNA), color=pRp_RNA)) + 
  geom_point() + ggtitle("Ribsonmal gene percentage") + ylab("log10(feature count)") + xlab("log10(umi count)") +
  geom_vline(aes(xintercept=log10(nCount_RNA_min)), color="red", linetype="longdash") +
  geom_vline(aes(xintercept=log10(nCount_RNA_max)), color="red", linetype="longdash") +
  geom_hline(aes(yintercept=log10(nFeature_RNA_min)), color="red", linetype="longdash") + 
  geom_hline(aes(yintercept=log10(nFeature_RNA_max)), color="red", linetype="longdash") + 
  facet_grid(tissue~treatment+sample_rep) + theme(aspect.ratio=1, legend.position="bottom") + 
  scale_size(guide=guide_legend(direction="vertical"))

options(repr.plot.width=15, repr.plot.height=10)
sc_1 + sc_2 + sc_3 + plot_layout(ncol=2) & theme(legend.position="bottom")
ggsave(sc_1, filename="result/plot/seurat/sc_mt.png", width=5, height=5)
ggsave(sc_2, filename="result/plot/seurat/sc_hg.png", width=5, height=5)
ggsave(sc_3, filename="result/plot/seurat/sc_rb.png", width=5, height=5)

so_qc <- subset(so_qc, subset=nFeature_RNA >= nFeature_RNA_min & nFeature_RNA <= nFeature_RNA_max & nCount_RNA >= nCount_RNA_min & nCount_RNA <= nCount_RNA_max)

Correlate percent Mt with OXPHOS nDNA model

Frederick et al., 2020 https://www.nature.com/articles/s41598-020-63448-z

go_bp <- msigdbr(species="mouse", category="C5", subcategory="BP")

# Positive regulation of glycolytic process (GO:0045821)
go_bp_gly <- go_bp[go_bp$gs_exact_source == "GO:0045821", ]$gene_symbol

# Positive regulation of oxidative phosphorylation (GO:1903862)
# Not found in msigdbr getting genes manually from http://www.informatics.jax.org/go/term/GO:1903862
go_bp_op <- c("Gadd45gip1", "Lexm", "Myc", "Myog", "Nupr1", "Vcp")
go_bp_op <- go_bp_op[go_bp_op %in% rownames(so_qc)]

options(warn=-1)
so_qc <- AddModuleScore(so_qc, features=list(go_bp_gly), assays="RNA", slot="data", ctrl=100, nbin=25, name="msGly_RNA")
so_qc <- AddModuleScore(so_qc, features=list(go_bp_op), assays="RNA", slot="data", ctrl=100, nbin=25, name="msOx_RNA")

sc_1 <- ggplot(so_qc@meta.data, aes(x=pMt_RNA, y=msGly_RNA1, color=log10(nFeature_RNA))) + 
    geom_point() + ggtitle("Glycolytic modele score vs %MT") + ylab("Module score Glycolysis") + xlab("%Mt") +
    geom_vline(aes(xintercept=pMt_RNA_max), color="red", linetype="longdash") +
    facet_grid(tissue~treatment+sample_rep) + theme(aspect.ratio=1, legend.position="bottom") + 
    scale_size(guide=guide_legend(direction="vertical"))

sc_2 <- ggplot(so_qc@meta.data, aes(x=pMt_RNA, y=msOx_RNA1, color=log10(nFeature_RNA))) + 
    geom_point() + ggtitle("OXPHOS module score vs %MT") + ylab("Module score OXPHOS") + xlab("%Mt") +
    geom_vline(aes(xintercept=pMt_RNA_max), color="red", linetype="longdash") +
    facet_grid(tissue~treatment+sample_rep) + theme(aspect.ratio=1, legend.position="bottom") + 
    scale_size(guide=guide_legend(direction="vertical"))

options(repr.plot.width=15, repr.plot.height=5)
ggsave(sc_1, filename="result/plot/seurat/sc_mt_glycolytic.png", width=5, height=5)
ggsave(sc_2, filename="result/plot/seurat/sc_mt_oxphos.png", width=5, height=5)
sc_1 + sc_2 + plot_layout(ncol=2) & theme(legend.position="bottom")

so_qc <- subset(so_qc, subset=pMt_RNA <= pMt_RNA_max)

Box plots

qc_vln_FUN <- function(data, y, fill, ylab="", ymin, ymax) {
  
  vln_plot_1 <- ggplot(data, aes(x=treatment, y={{y}}, color={{fill}})) + 
    geom_jitter(alpha=0.2, shape=16, color="gray") + 
    geom_boxplot(alpha=1.0) + xlab("") + 
    ylim(ymin, ymax) +
    scale_color_manual(values=color$tissue) + 
    ggtitle("GEM") + 
    facet_wrap(~tissue, scales="free_x") + ylab(ylab) +
    theme(
      plot.title=element_text(size=12, face="bold", margin=margin(t=0, r=0, b=5, l=0)), 
      axis.text.x=element_text(angle=45, vjust=1, hjust=1), 
      strip.text=element_blank()
      )
  
  vln_plot_2 <- ggplot(data[data$qc_class == "pass", ], aes(x=treatment, y={{y}}, color={{fill}})) +
    geom_jitter(alpha=0.2, shape=16, color="gray") +
    geom_boxplot(alpha=1.0) + xlab("") + 
    ylim(ymin, ymax) +
    scale_color_manual(values=color$tissue) +
    ggtitle("Filtered") +
    facet_wrap(~tissue, scales="free_x") + ylab(ylab) + 
    theme(
      plot.title=element_text(size=12, face="bold", margin=margin(t=0, r=0, b=5, l=0)), 
      axis.text.x=element_text(angle=45, vjust=1, hjust=1), 
      strip.text=element_blank()
      )

  return(list(vln_plot_1, vln_plot_2))
}

data <- dplyr::filter(so_raw@meta.data, cellranger_class == "Cell")

qc_vln_1 <- qc_vln_FUN(data, nCount_RNA, tissue, "UMI [count]",  ymin=0, ymax=max(so_raw$nCount_RNA)+100)
qc_vln_2 <- qc_vln_FUN(data, nFeature_RNA, tissue, "Feature [count]", ymin=0, ymax=max(so_raw$nFeature_RNA)+100)
qc_vln_3 <- qc_vln_FUN(data, pMt_RNA, tissue, "Mt [%]", ymin=0, ymax=100)
qc_vln_4 <- qc_vln_FUN(data, pHb_RNA, tissue, "Hb [%]", ymin=0, ymax=100)
qc_vln_5 <- qc_vln_FUN(data, pRp_RNA, tissue, "Rbl [%]", ymin=0, ymax=100)

options(repr.plot.width=15, repr.plot.height=10)
qc_vln_1[[1]] + qc_vln_1[[2]] + qc_vln_2[[1]]  + qc_vln_2[[2]] + 
qc_vln_3[[1]] + qc_vln_3[[2]] + qc_vln_4[[1]]  + qc_vln_4[[2]] + 
qc_vln_5[[1]] + qc_vln_5[[2]] + plot_spacer() + plot_spacer() + plot_layout(guides="collect", ncol=4) & theme(legend.position="bottom")

ggsave(qc_vln_1[[1]] + qc_vln_1[[2]] + plot_layout(guides="collect", ncol=2) & theme(legend.position="bottom"), filename="result/plot/seurat/qc_vln_1.png", width=5, height=2.5)
ggsave(qc_vln_2[[1]] + qc_vln_2[[2]] + plot_layout(guides="collect", ncol=2) & theme(legend.position="bottom"), filename="result/plot/seurat/qc_vln_2.png", width=5, height=2.5)
ggsave(qc_vln_3[[1]] + qc_vln_3[[2]] + plot_layout(guides="collect", ncol=2) & theme(legend.position="bottom"), filename="result/plot/seurat/qc_vln_3.png", width=5, height=2.5)
ggsave(qc_vln_4[[1]] + qc_vln_4[[2]] + plot_layout(guides="collect", ncol=2) & theme(legend.position="bottom"), filename="result/plot/seurat/qc_vln_4.png", width=5, height=2.5)
ggsave(qc_vln_5[[1]] + qc_vln_5[[2]] + plot_layout(guides="collect", ncol=2) & theme(legend.position="bottom"), filename="result/plot/seurat/qc_vln_5.png", width=5, height=2.5)

SingleR annotation

SingleR identifies marker genes from the reference and uses them to compute assignment scores (based on the Spearman correlation across markers) for each cell in the test dataset against each label in the reference. The label with the highest score is the assigned to the test cell, possibly with further fine-tuning to resolve closely related labels.

first.labels: Labels before fine-tuning
labels: Labels after fine-tuning
pruning: labels after pruning

ref <- readRDS("data/haemosphere/se_haemosphere.rds")

# Seurat object to SingleCellExperiment
sce <- SingleCellExperiment(list(counts=so_qc@assays$RNA@counts))
    
# Predict labels
label_main <- SingleR::SingleR(test=sce, ref=ref, labels=ref$label.main, assay.type.test="counts", de.method="classic")
label_fine <- SingleR::SingleR(test=sce, ref=ref, labels=ref$label.fine, assay.type.test="counts", de.method="classic")

# Save predicted labels 
saveRDS(label_main, "result/label_main.rds")
saveRDS(label_fine, "result/label_fine.rds")

# Add labels to Seurat object
label_main_meta <- as.data.frame(label_main) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(main_labels=pruned.labels, main_delta_score=tuning.scores.first)
label_fine_meta <- as.data.frame(label_fine) %>% dplyr::select(pruned.labels, tuning.scores.first) %>% dplyr::rename(fine_labels=pruned.labels, fine_delta_score=tuning.scores.first)

so_qc <- AddMetaData(so_qc, cbind(label_main_meta, label_fine_meta))

# Set factor level for labels 
so_qc$main_labels <- factor(so_qc$main_labels, levels=names(color$main_labels_haemosphere))
so_qc$fine_labels <- factor(so_qc$fine_labels, levels=names(color$fine_labels_haemosphere))

so_qc$main_labels <- factor(so_qc$main_labels, levels=names(color$main_labels_haemosphere))
so_qc$fine_labels <- factor(so_qc$fine_labels, levels=names(color$fine_labels_haemosphere))

Ery score

so_qc$ery_score=label_main$scores[, "Ery"]

so_qc$main_labels <- factor(so_qc$main_labels, levels=names(color$main_labels_haemosphere))

boxplot_ery_score <- ggplot(so_qc@meta.data, aes(x=ery_score, color=main_labels)) + 
    geom_boxplot() + 
    scale_color_manual(values=color$main_labels_haemosphere) + 
    guides(colour=guide_legend(reverse=TRUE)) + 
    facet_grid(~tissue, scales="free_y")

count_module <- function(so, pattern, assay="RNA", slot="counts", stat="sum") {
    
    genes <- rownames(so)[grep(pattern, rownames(so))]
    cnt <- GetAssayData(so, assay=assay, slot=slot)
    cnt <- as.matrix(cnt)
    cnt <- cnt[rownames(cnt) %in% genes, ]
    
    if(stat=="sum"){cnt_stat <- colSums(cnt)}

    return(cnt_stat)
}
so_qc$cHb_RNA <- count_module(so_qc, pattern="Hba-|Hbb-|Hbq1b|Hbq1a")

boxplot_ery_count_module <- ggplot(so_qc@meta.data, aes(x=log10(cHb_RNA+1), color=main_labels)) + 
    geom_boxplot() + 
    scale_color_manual(values=color$main_labels_haemosphere) + 
    guides(colour=guide_legend(reverse=TRUE)) +
    facet_grid(~tissue, scales="free_y")

boxplot_ery_count_feature <- ggplot(so_qc@meta.data, aes(x=nFeature_RNA, y=log10(cHb_RNA), color=main_labels)) + 
    geom_point() + 
    geom_hline(aes(yintercept=log10(10)), color="red", linetype="longdash") +
    scale_color_manual(values=color$main_labels_haemosphere) + 
    guides(colour=guide_legend(reverse=TRUE)) +
    facet_grid(~tissue, scales="free_y")

boxplot_ery_count_umi <- ggplot(so_qc@meta.data, aes(x=nCount_RNA, y=log10(cHb_RNA), color=main_labels)) + 
    geom_point() + 
    geom_hline(aes(yintercept=log10(10)), color="red", linetype="longdash") +
    scale_color_manual(values=color$main_labels_haemosphere) + 
    guides(colour=guide_legend(reverse=TRUE)) +
    facet_grid(~tissue, scales="free_y")

options(repr.plot.width=10, repr.plot.height=10, warn=-1)
boxplot_ery_score + boxplot_ery_count_module + boxplot_ery_count_feature + boxplot_ery_count_umi + plot_layout(guides="collect", ncol=2) & theme(legend.position="bottom")

main_labels_split <- split(label_main, f=so_qc$sample_group)
main_labels_score_hm <- list()
for(i in names(main_labels_split)) {main_labels_score_hm[[i]] <- singler_score_hm(main_labels_split[[i]], color=color$main_labels_haemosphere, main=i)}

options(repr.plot.width=20, repr.plot.height=20)
main_labels_score_hm_grid <- gridExtra::arrangeGrob(grobs=main_labels_score_hm, ncol=1) 
ggsave(main_labels_score_hm_grid, filename="result/plot/seurat/main_labels_score_hm_grid.png", width=20, height=20)
grid::grid.draw(main_labels_score_hm_grid)

so_qc$ery_qc <- ifelse(!(so_qc$tissue == "Myeloid" & so_qc$main_labels == "Ery"), "pass", "fail")
so_qc <- subset(so_qc, subset = ery_qc == "pass")

Compute module scores

Histocompatibility Class I

H2-K1 - Class I histocompatibility antigen, kappa-B alpha chain
H2-D1 - Class I histocompatibility antigen, D-B alpha chain
H2-L1 - not found

Histocompatibility Class II

H2-Aa - Class II histocompatibility antigen, A-B alpha chain
H2-Ab1 - Class II histocompatibility antigen, A beta chain

H2-Eb1 - Class II histocompatibility antigen, I-E beta chain
H2-Eb2 - Class II histocompatibility antigen, E beta chain

Info: The antigens presented by class I peptides are derived from cytosolic proteins while calss II derived from extracellular proteins.

mhcI_genes <- c("H2-K1", "H2-D1")
mhcII_genes <- c("H2-Aa", "H2-Ab1", "H2-Eb1", "H2-Eb2")

so_qc <- AddModuleScore(so_qc, features=list(mhcI_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msMHCI_RNA")
so_qc <- AddModuleScore(so_qc, features=list(mhcII_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msMHCII_RNA")

Surface marker

Cd4 - T-cell surface glycoprotein CD4, Co-receptor for MHC class II, active on T-helper cells and triggers differentiation of monocytes into functional mature macrophages.
Cd8 - T-cell surface glycoprotein CD8, Co-receptor for MHC class I, CTL activattion and thymic selection of Cd8+ T cells.

cd4_genes <- c("Cd4") 
cd8_genes <- c("Cd8a", "Cd8b1")

so_qc <- AddModuleScore(so_qc, features=list(cd4_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msCd4_RNA")
so_qc <- AddModuleScore(so_qc, features=list(cd8_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msCd8_RNA")

T-cell TCR-CD3 complex

Trac - T cell receptor alpha constant
Trbc1 - T cell receptor beta constant 1
Trbc2 - T cell receptor beta constant 2
Trdc - T cell receptor delta constant
Trgc1 - T cell receptor gamma constant 1
Trgc2 - T cell receptor gamma constant 2

Cd247 - T-cell surface glycoprotein Cd3 zeta chain (Cd3z)
Cd3g - T-cell surface glycoprotein Cd3 gamma chain
Cd3e - T-cell surface glycoprotein Cd3 epsilon chain
Cd3d - T-cell surface glycoprotein Cd3 delta chain

tcr_genes <- c("Trac", "Trbc1", "Trbc2", "Trdc", "Trgc1", "Trgc2")
tcr_cd3_genes <- c("Cd247", "Cd3g", "Cd3e", "Cd3d")

so_qc <- AddModuleScore(so_qc, features=list(tcr_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msTcr_RNA")
so_qc <- AddModuleScore(so_qc, features=list(tcr_cd3_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msTcr_cd3_RNA")

B-cell Immonogloblin complex

Naive B-cells produce the following Ig classes:
Ighm - Immunoglobulin heavy constant mu (naive B-cells)
Ighd - Immunoglobulin heavy constant delta (naive B-cells)

Through isotope switching the following Ig classes can be produced:
Ighg1 - Immunoglobulin heavy constant gamma (Mouse with Igh1b have Igg2c isotope instead of Igg2a)
Ighg2a - Immunoglobulin heavy constant gamma (NA)
Ighg2b - Immunoglobulin heavy constant gamma
Ighg2c - Immunoglobulin heavy constant gamma
Ighg3 - Immunoglobulin heavy constant gamma
Igha - Immunoglobulin heavy constant alpha
Ighe - Immunoglobulin heavy constant epsilon (NA)

Igkc - Immunoglobulin kappa constant (light chain)
Iglc - Immunoglobulin lambda constant (light chain)

Ighm_genes <- c("Ighm")
Ighd_genes <- c("Ighd")
Ighg_genes <- c("Ighg1", "Ighg2a", "Ighg2b", "Ighg2c", "Ighg3")
Igha_genes <- c("Igha")
Ighe_genes <- c("Ighe") 
Igkc_genes <- c("Igkc")
Iglc_genes <- c("Iglc1", "Iglc2", "Iglc3")

so_qc <- AddModuleScore(so_qc, features=list(Ighm_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msIghm_RNA")
so_qc <- AddModuleScore(so_qc, features=list(Ighd_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msIghd_RNA")
so_qc <- AddModuleScore(so_qc, features=list(Ighg_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msIghg_RNA")
so_qc <- AddModuleScore(so_qc, features=list(Igha_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msIgha_RNA")
#so_qc <- AddModuleScore(so_qc, features=list(Ighe_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msIghe_RNA")
so_qc <- AddModuleScore(so_qc, features=list(Igkc_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msIgkc_RNA")
so_qc <- AddModuleScore(so_qc, features=list(Iglc_genes), assays="RNA", slot="data", ctrl=100, nbin=25, name="msIglc_RNA")

Compute cell cycle score

cc_genes_seurat_s <- cc.genes.updated.2019$s.genes
cc_genes_seurat_g2m <- cc.genes.updated.2019$g2m.genes

# Translate human to mouse 
library(biomaRt) 

# Set ssl for biomart  
httr::set_config(httr::config(ssl_verifypeer=FALSE))

ensembl_hgnc <- useMart("ensembl", dataset="hsapiens_gene_ensembl") 
ensembl_mm <- useMart("ensembl", dataset="mmusculus_gene_ensembl")

cc_genes_seurat_s <- getLDS(attributes=c("hgnc_symbol"), filters="hgnc_symbol", values=cc_genes_seurat_s, mart=ensembl_hgnc, attributesL=c("mgi_symbol"), martL=ensembl_mm, uniqueRows=T)
cc_genes_seurat_s <- cc_genes_seurat_s[, 2]

cc_genes_seurat_g2m <- getLDS(attributes=c("hgnc_symbol"), filters="hgnc_symbol", values=cc_genes_seurat_g2m, mart=ensembl_hgnc, attributesL=c("mgi_symbol"), martL=ensembl_mm, uniqueRows=T)
cc_genes_seurat_g2m <- cc_genes_seurat_g2m[, 2]

Ensembl site unresponsive, trying asia mirror

so_qc <- CellCycleScoring(so_qc, s.features=cc_genes_seurat_s, g2m.features=cc_genes_seurat_g2m, set.ident=FALSE)

colnames(so_qc@meta.data) <- gsub("Phase", "cc_phase_class", colnames(so_qc@meta.data))
colnames(so_qc@meta.data) <- gsub("S.Score", "msS_RNA", colnames(so_qc@meta.data))
colnames(so_qc@meta.data) <- gsub("G2M.Score", "msG2M_RNA", colnames(so_qc@meta.data)) 
    
so_qc$msCC_diff_RNA <- so_qc$msS_RNA - so_qc$msG2M_RNA
so_qc$cc_phase_class <- factor(so_qc$cc_phase_class, levels=names(color$cc_phase_class))

saveRDS(so_qc, "data/object/so_qc.rds")

SCTransform and dimensional reduction (sample_group)

options(warn=-1)
so_qc_sample_group <- SplitObject(so_qc, split.by="sample_group")
so_qc_sample_group <- lapply(so_qc_sample_group, function(so) {
    
    so <- SCTransform(so, assay="RNA", verbose=FALSE)
    so <- RunPCA(so, npcs=100, verbose=FALSE)
    so <- FindNeighbors(so, dims=1:100, k.param=10, verbose=FALSE)
    so <- FindClusters(so, verbose=FALSE)
    so <- RunUMAP(so, dims=1:100, verbose=FALSE)
    
    return(so)})

saveRDS(so_qc_sample_group, "data/object/so_qc_sample_group.rds")

options(repr.plot.width=20, repr.plot.height=20)
dim_plot_sample_group <- lapply(so_qc_sample_group, dim_plot, cluster="SCT_snn_res.0.8")

ggsave(dim_plot_sample_group[[1]], filename="result/plot/seurat/dim_plot_sample_group_m_cpg.png", width=20, height=20)
ggsave(dim_plot_sample_group[[2]], filename="result/plot/seurat/dim_plot_sample_group_p_cpg.png", width=20, height=20)
ggsave(dim_plot_sample_group[[3]], filename="result/plot/seurat/dim_plot_sample_group_m_nacl.png", width=20, height=20)
ggsave(dim_plot_sample_group[[4]], filename="result/plot/seurat/dim_plot_sample_group_p_nacl.png", width=20, height=20)

options(repr.plot.width=20, repr.plot.height=10)

so_qc_sample_group_pc_annotation <- lapply(so_qc_sample_group, pc_annotation, dims=100, ident="sample_group")
so_qc_sample_group_pc_annotation_plot <- lapply(so_qc_sample_group_pc_annotation, pc_annotation_plot)

so_qc_sample_group_pc_annotation_plot <- 
    so_qc_sample_group_pc_annotation_plot[[1]] + 
    so_qc_sample_group_pc_annotation_plot[[2]] + 
    so_qc_sample_group_pc_annotation_plot[[3]] + 
    so_qc_sample_group_pc_annotation_plot[[4]] + 
    plot_layout(ncol=2, guides="collect")

so_qc_sample_group_pc_annotation_plot
ggsave(so_qc_sample_group_pc_annotation_plot, filename="result/plot/seurat/so_qc_sample_group_pc_loadings_plot.png", width=10, height=5)

Identify outlier cluster

Sample replicates

options(repr.plot.width=10, repr.plot.height=10, warn=-1)
cluster_sample_rep <- function(so) {
    
    box_plot <- ggplot(so@meta.data, aes(x=SCT_snn_res.0.8, fill=tissue, alpha=sample_rep)) + 
        geom_bar(stat="count", position="fill") + 
        scale_fill_manual(values=color$tissue) +
        scale_alpha_manual(values=c(1, 0.7)) + 
        ggtitle(paste0(so$tissue[1], so$treatment[1])) + xlab("Cluster") + ylab("Cell frequency") + 
        facet_grid(~tissue)
    
    return(box_plot)
}

cluster_sample_rep_plot <- lapply(so_qc_sample_group, cluster_sample_rep)
cluster_sample_rep_plot <- 
    cluster_sample_rep_plot[[1]] + 
    cluster_sample_rep_plot[[2]] + 
    cluster_sample_rep_plot[[3]] + 
    cluster_sample_rep_plot[[4]] + 
    plot_layout(ncol=2, guides="collect") & theme(legend.position="bottom")

cluster_sample_rep_plot

Myeloid CpG

# Myeloid_CpG
so_qc_sample_group_m_cpg <- so_qc_sample_group[[1]]

# Cluster 16: B-cell 
deg_cluster_16 <- FindMarkers(so_qc_sample_group_m_cpg, ident.1=16, only.pos=FALSE, group.by="SCT_snn_res.0.8")

options(repr.plot.width=20, repr.plot.height=5)
deg_vp_sample_group_m_cpg <- 
    deg_volcano_plot(deg_cluster_16, title="Myeloid (CpG) - Cluster 16") + 
    plot_spacer() + plot_spacer() + plot_spacer() +
    plot_layout(ncol=4)
deg_vp_sample_group_m_cpg
ggsave(deg_vp_sample_group_m_cpg, filename="result/plot/seurat/deg_vp_sample_group_m_cpg.png", width=20, height=5)

Progenitor CpG

# Progenitor CpG
so_qc_sample_group_p_cpg <- so_qc_sample_group[[2]]

# Cluster 14: B-cell 
deg_cluster_14 <- FindMarkers(so_qc_sample_group_p_cpg, ident.1=14, only.pos=FALSE, group.by="SCT_snn_res.0.8")

options(repr.plot.width=20, repr.plot.height=5)
deg_vp_sample_group_p_cpg <- 
    deg_volcano_plot(deg_cluster_14, title="Progenitor (CpG) - Cluster 14") +
    plot_spacer() + plot_spacer() + plot_spacer() +
    plot_layout(ncol=4)
deg_vp_sample_group_p_cpg
ggsave(deg_vp_sample_group_p_cpg, filename="result/plot/seurat/deg_vp_sample_group_p_cpg.png", width=20, height=5)

Myeloid NaCl

# Myeloid NaCl
so_qc_sample_group_m_nacl <- so_qc_sample_group[[3]]

# Cluster 17: B-cells 
deg_cluster_17 <- FindMarkers(so_qc_sample_group_m_nacl, ident.1=17, only.pos=FALSE, group.by="SCT_snn_res.0.8")

# Cluster 10, 12, 13, 15, 16: Unknown
deg_cluster_10 <- FindMarkers(so_qc_sample_group_m_nacl, ident.1=10, only.pos=FALSE, group.by="SCT_snn_res.0.8")
deg_cluster_12 <- FindMarkers(so_qc_sample_group_m_nacl, ident.1=12, only.pos=FALSE, group.by="SCT_snn_res.0.8")
deg_cluster_13 <- FindMarkers(so_qc_sample_group_m_nacl, ident.1=13, only.pos=FALSE, group.by="SCT_snn_res.0.8")
deg_cluster_15 <- FindMarkers(so_qc_sample_group_m_nacl, ident.1=15, only.pos=FALSE, group.by="SCT_snn_res.0.8")
deg_cluster_16 <- FindMarkers(so_qc_sample_group_m_nacl, ident.1=16, only.pos=FALSE, group.by="SCT_snn_res.0.8")

options(repr.plot.width=20, repr.plot.height=10)
deg_vp_sample_group_m_nacl <- 
    deg_volcano_plot(deg_cluster_10, title="Myeloid (NaCl) - Cluster 10") +
    deg_volcano_plot(deg_cluster_12, title="Myeloid (NaCl) - Cluster 12") +
    deg_volcano_plot(deg_cluster_13, title="Myeloid (NaCl) - Cluster 13") +
    deg_volcano_plot(deg_cluster_15, title="Myeloid (NaCl) - Cluster 15") +
    deg_volcano_plot(deg_cluster_16, title="Myeloid (NaCl) - Cluster 16") +
    deg_volcano_plot(deg_cluster_17, title="Myeloid (NaCl) - Cluster 17") +
    plot_spacer() + plot_spacer() +    
    plot_layout(ncol=4)
deg_vp_sample_group_m_nacl
ggsave(deg_vp_sample_group_m_nacl, filename="result/plot/seurat/deg_vp_sample_group_m_nacl.png", width=20, height=10)

Progenitor NaCl

# Myeloid NaCl
so_qc_sample_group_p_nacl <- so_qc_sample_group[[4]]

# Cluster 11: T-cell  
deg_cluster_11 <- FindMarkers(so_qc_sample_group_p_nacl, ident.1=11, only.pos=FALSE, group.by="SCT_snn_res.0.8")

# Cluster 15: B-cell  
deg_cluster_15 <- FindMarkers(so_qc_sample_group_p_nacl, ident.1=15, only.pos=FALSE, group.by="SCT_snn_res.0.8")

options(repr.plot.width=20, repr.plot.height=5)
deg_vp_sample_group_p_nacl <- 
    deg_volcano_plot(deg_cluster_11, title="Progenitor (NaCl) - Cluster 11") +
    deg_volcano_plot(deg_cluster_15, title="Progenitor (NaCl) - Cluster 15") +
    plot_spacer() + plot_spacer() + 
    plot_layout(ncol=4)
deg_vp_sample_group_p_nacl
ggsave(deg_vp_sample_group_p_nacl, filename="result/plot/seurat/deg_vp_sample_group_p_nacl.png", width=20, height=5)

Remove cluster from so_qc (sample group)

so_qc <- merge(so_qc_sample_group[[1]], so_qc_sample_group[2:4])

so_qc <- subset(so_qc, subset=!(
    (tissue == "Myeloid" & treatment == "CpG" &  SCT_snn_res.0.8 %in% c(16)) | 
    (tissue == "Progenitor" & treatment == "CpG" & SCT_snn_res.0.8 %in% c(14)) | 
    (tissue == "Myeloid" & treatment == "NaCl" & SCT_snn_res.0.8 %in% c(15, 13, 17)) | 
    (tissue == "Progenitor" & treatment == "NaCl" & SCT_snn_res.0.8 %in% c(11, 15))
)
      )

SCTransform and dimensional reduction (treatment)

options(warn=-1)
so_qc_treatment <- SplitObject(so_qc, split.by="treatment")
so_qc_treatment <- lapply(so_qc_treatment, function(so) {
    
    so <- SCTransform(so, assay="RNA", verbose=FALSE)
    so <- RunPCA(so, npcs=100, verbose=FALSE)
    so <- FindNeighbors(so, dims=1:100, verbose=FALSE)
    so <- FindClusters(so, verbose=FALSE)
    so <- RunUMAP(so, dims=1:100, verbose=FALSE)
    
    return(so)})

saveRDS(so_qc_treatment, "data/object/so_qc_treatment.rds")

options(repr.plot.width=20, repr.plot.height=20)
dim_plot_treatment <- lapply(so_qc_treatment, dim_plot, cluster="SCT_snn_res.0.8")

ggsave(dim_plot_treatment[[1]], filename="result/plot/seurat/dim_plot_treatment_cpg.png", width=20, height=20)
ggsave(dim_plot_treatment[[2]], filename="result/plot/seurat/dim_plot_treatment_nacl.png", width=20, height=20)

options(repr.plot.width=20, repr.plot.height=5)

so_qc_treatment_pc_annotation <- lapply(so_qc_treatment, pc_annotation, dims=100, ident="treatment")
so_qc_treatment_pc_annotation_plot <- lapply(so_qc_treatment_pc_annotation, pc_annotation_plot)

so_qc_treatment_pc_annotation_plot <- 
    so_qc_treatment_pc_annotation_plot[[1]] + 
    so_qc_treatment_pc_annotation_plot[[2]] + 
    plot_layout(ncol=2, guides="collect")

so_qc_treatment_pc_annotation_plot
ggsave(so_qc_treatment_pc_annotation_plot, filename="result/plot/seurat/so_qc_treatment_pca_loadings_plot.png", width=10, height=5)

options(repr.plot.width=20, repr.plot.height=5, warn=-1)
cluster_sample_rep <- function(so) {
    
    box_plot <- ggplot(so@meta.data, aes(x=SCT_snn_res.0.8, fill=tissue, alpha=sample_rep)) + 
        geom_bar(stat="count", position="fill") + 
        scale_fill_manual(values=color$tissue) +
        scale_alpha_manual(values=c(1, 0.7)) + 
        ggtitle(so$treatment[1]) + xlab("Cluster") + ylab("Cell frequency") + 
        facet_grid(~tissue)
    
    return(box_plot)
}

cluster_sample_rep_plot <- lapply(so_qc_treatment, cluster_sample_rep)
cluster_sample_rep_plot <- 
    cluster_sample_rep_plot[[1]] + 
    cluster_sample_rep_plot[[2]] + 
    plot_layout(ncol=2, guides="collect") & theme(legend.position="bottom")
cluster_sample_rep_plot

options(repr.plot.width=10, repr.plot.height=5)

cluster_tissue_1 <- ggplot(so_qc_treatment[[1]]@meta.data, aes(x=SCT_snn_res.0.8, fill=tissue)) + 
      geom_bar(stat="count", position="fill") + 
      scale_fill_manual(values=color$tissue) +
      ggtitle("CpG sample") + xlab("Cluster") + ylab("Cell frequency") + 
      theme(legend.position="bottom")

cluster_tissue_2 <- ggplot(so_qc_treatment[[2]]@meta.data, aes(x=SCT_snn_res.0.8, fill=tissue)) + 
      geom_bar(stat="count", position="fill") + 
      scale_fill_manual(values=color$tissue) +
      ggtitle("NaCl sample") + xlab("Cluster") + ylab("Cell frequency") + 
      theme(legend.position="bottom")

cluster_tissue <- cluster_tissue_1 + cluster_tissue_2
cluster_tissue
ggsave(cluster_tissue, filename="result/plot/seurat/cluster_tissue_treatment.png", width=10, height=2.5)

SCTransform + Cell cycle Regression and dimensional reduction (treatment)

options(warn=-1)
so_qc_treatment_reg <- SplitObject(so_qc, split.by="treatment")
so_qc_treatment_reg <- lapply(so_qc_treatment_reg, function(so) {
    
    so <- SCTransform(so, assay="RNA", vars.to.regress=c("msCC_diff_RNA"), verbose=FALSE)
    so <- RunPCA(so, npcs=100, verbose=FALSE)
    so <- FindNeighbors(so, dims=1:100, verbose=FALSE)
    so <- FindClusters(so, verbose=FALSE)
    so <- RunUMAP(so, dims=1:100, verbose=FALSE)
    
    return(so)})

saveRDS(so_qc_treatment_reg, "data/object/so_qc_treatment_reg.rds")

options(repr.plot.width=20, repr.plot.height=20)
dim_plot_treatment_reg <- lapply(so_qc_treatment_reg, dim_plot, cluster="SCT_snn_res.0.8")

ggsave(dim_plot_treatment_reg[[1]], filename="result/plot/seurat/dim_plot_treatment_reg_cpg.png", width=20, height=20)
ggsave(dim_plot_treatment_reg[[2]], filename="result/plot/seurat/dim_plot_treatment_reg_nacl.png", width=20, height=20)

options(repr.plot.width=20, repr.plot.height=5)

so_qc_treatment_reg_pc_annotation <- lapply(so_qc_treatment_reg, pc_annotation, dims=100, ident="treatment")
so_qc_treatment_reg_pc_annotation_plot <- lapply(so_qc_treatment_reg_pc_annotation, pc_annotation_plot)

so_qc_treatment_reg_pc_annotation_plot <- 
    so_qc_treatment_reg_pc_annotation_plot[[1]] + 
    so_qc_treatment_reg_pc_annotation_plot[[2]] + 
    plot_layout(ncol=2, guides="collect")

so_qc_treatment_reg_pc_annotation_plot
ggsave(so_qc_treatment_reg_pc_annotation_plot, filename="result/plot/seurat/so_qc_treatment_reg_pca_loadings_plot.png", width=10, height=5)

options(repr.plot.width=20, repr.plot.height=5, warn=-1)
cluster_sample_rep <- function(so) {
    
    box_plot <- ggplot(so@meta.data, aes(x=SCT_snn_res.0.8, fill=tissue, alpha=sample_rep)) + 
        geom_bar(stat="count", position="fill") + 
        scale_fill_manual(values=color$tissue) +
        scale_alpha_manual(values=c(1, 0.7)) + 
        ggtitle(so$treatment[1]) + xlab("Cluster") + ylab("Cell frequency") + 
        facet_grid(~tissue)
    
    return(box_plot)
}

cluster_sample_rep_plot <- lapply(so_qc_treatment_reg, cluster_sample_rep)
cluster_sample_rep_plot <- 
    cluster_sample_rep_plot[[1]] + 
    cluster_sample_rep_plot[[2]] + 
    plot_layout(ncol=2, guides="collect") & theme(legend.position="bottom")
cluster_sample_rep_plot

options(repr.plot.width=10, repr.plot.height=5)
cluster_tissue_1 <- ggplot(so_qc_treatment_reg[[1]]@meta.data, aes(x=SCT_snn_res.0.8, fill=tissue)) + 
      geom_bar(stat="count", position="fill") + 
      scale_fill_manual(values=color$tissue) +
      ggtitle("CpG sample") + xlab("Cluster") + ylab("Cell frequency") + 
      theme(legend.position="bottom")

cluster_tissue_2 <- ggplot(so_qc_treatment_reg[[2]]@meta.data, aes(x=SCT_snn_res.0.8, fill=tissue)) + 
      geom_bar(stat="count", position="fill") + 
      scale_fill_manual(values=color$tissue) +
      ggtitle("NaCl sample") + xlab("Cluster") + ylab("Cell frequency") + 
      theme(legend.position="bottom")

cluster_tissue <- cluster_tissue_1 + cluster_tissue_2
cluster_tissue
ggsave(cluster_tissue, filename="result/plot/seurat/cluster_tissue_treatment_reg.png", width=10, height=5)

Save output

# saveRDS(so_qc_treatment_reg, so_qc_file)

Session info

sessionInfo()

R version 4.1.0 (2021-05-18)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /home/fdeckert/bin/miniconda3/envs/r.4.1.0-FD20200109SPLENO/lib/libopenblasp-r0.3.15.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
 [1] grid      parallel  stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] ggrepel_0.9.1               biomaRt_2.48.2             
 [3] viridis_0.6.1               viridisLite_0.4.0          
 [5] RColorBrewer_1.1-2          patchwork_1.1.1            
 [7] gridExtra_2.3               ComplexHeatmap_2.8.0       
 [9] ggplot2_3.3.3               msigdbr_7.4.1              
[11] Matrix_1.3-4                reshape2_1.4.4             
[13] dplyr_1.0.6                 celldex_1.2.0              
[15] SingleR_1.6.1               SingleCellExperiment_1.14.1
[17] SummarizedExperiment_1.22.0 Biobase_2.52.0             
[19] GenomicRanges_1.44.0        GenomeInfoDb_1.28.0        
[21] IRanges_2.26.0              S4Vectors_0.30.0           
[23] BiocGenerics_0.38.0         MatrixGenerics_1.4.0       
[25] matrixStats_0.59.0          SeuratObject_4.0.1         
[27] Seurat_4.0.2               

loaded via a namespace (and not attached):
  [1] utf8_1.2.1                    reticulate_1.20              
  [3] tidyselect_1.1.1              RSQLite_2.2.5                
  [5] AnnotationDbi_1.54.0          htmlwidgets_1.5.3            
  [7] BiocParallel_1.26.0           Rtsne_0.15                   
  [9] munsell_0.5.0                 ScaledMatrix_1.0.0           
 [11] codetools_0.2-18              ica_1.0-2                    
 [13] pbdZMQ_0.3-5                  future_1.21.0                
 [15] miniUI_0.1.1.1                withr_2.4.2                  
 [17] colorspace_2.0-1              filelock_1.0.2               
 [19] uuid_0.1-4                    ROCR_1.0-11                  
 [21] tensor_1.5                    listenv_0.8.0                
 [23] labeling_0.4.2                repr_1.1.3                   
 [25] GenomeInfoDbData_1.2.6        polyclip_1.10-0              
 [27] farver_2.1.0                  bit64_4.0.5                  
 [29] parallelly_1.25.0             vctrs_0.3.8                  
 [31] generics_0.1.0                BiocFileCache_2.0.0          
 [33] R6_2.5.0                      doParallel_1.0.16            
 [35] clue_0.3-59                   rsvd_1.0.5                   
 [37] bitops_1.0-7                  spatstat.utils_2.1-0         
 [39] cachem_1.0.5                  DelayedArray_0.18.0          
 [41] assertthat_0.2.1              promises_1.2.0.1             
 [43] scales_1.1.1                  gtable_0.3.0                 
 [45] beachmat_2.8.0                Cairo_1.5-12.2               
 [47] globals_0.14.0                goftest_1.2-2                
 [49] rlang_0.4.11                  GlobalOptions_0.1.2          
 [51] splines_4.1.0                 lazyeval_0.2.2               
 [53] spatstat.geom_2.1-0           BiocManager_1.30.15          
 [55] yaml_2.2.1                    abind_1.4-5                  
 [57] httpuv_1.6.1                  tools_4.1.0                  
 [59] ellipsis_0.3.2                spatstat.core_2.1-2          
 [61] ggridges_0.5.3                Rcpp_1.0.6                   
 [63] plyr_1.8.6                    progress_1.2.2               
 [65] base64enc_0.1-3               sparseMatrixStats_1.4.0      
 [67] zlibbioc_1.38.0               purrr_0.3.4                  
 [69] RCurl_1.98-1.3                prettyunits_1.1.1            
 [71] rpart_4.1-15                  deldir_0.2-10                
 [73] pbapply_1.4-3                 GetoptLong_1.0.5             
 [75] cowplot_1.1.1                 zoo_1.8-9                    
 [77] cluster_2.1.2                 magrittr_2.0.1               
 [79] RSpectra_0.16-0               magick_2.7.2                 
 [81] data.table_1.14.0             scattermore_0.7              
 [83] circlize_0.4.13               lmtest_0.9-38                
 [85] RANN_2.6.1                    fitdistrplus_1.1-5           
 [87] hms_1.1.0                     mime_0.10                    
 [89] evaluate_0.14                 xtable_1.8-4                 
 [91] XML_3.99-0.6                  shape_1.4.6                  
 [93] compiler_4.1.0                tibble_3.1.2                 
 [95] KernSmooth_2.23-20            crayon_1.4.1                 
 [97] htmltools_0.5.1.1             mgcv_1.8-36                  
 [99] later_1.2.0                   tidyr_1.1.3                  
[101] DBI_1.1.1                     ExperimentHub_2.0.0          
[103] dbplyr_2.1.1                  MASS_7.3-54                  
[105] rappdirs_0.3.3                babelgene_21.4               
[107] igraph_1.2.6                  pkgconfig_2.0.3              
[109] IRdisplay_1.0                 plotly_4.9.3                 
[111] spatstat.sparse_2.0-0         xml2_1.3.2                   
[113] foreach_1.5.1                 XVector_0.32.0               
[115] stringr_1.4.0                 digest_0.6.27                
[117] sctransform_0.3.2             RcppAnnoy_0.0.18             
[119] spatstat.data_2.1-0           Biostrings_2.60.0            
[121] leiden_0.3.8                  uwot_0.1.10                  
[123] DelayedMatrixStats_1.14.0     curl_4.3.1                   
[125] shiny_1.6.0                   rjson_0.2.20                 
[127] lifecycle_1.0.0               nlme_3.1-152                 
[129] jsonlite_1.7.2                BiocNeighbors_1.10.0         
[131] limma_3.48.1                  fansi_0.5.0                  
[133] pillar_1.6.1                  lattice_0.20-44              
[135] KEGGREST_1.32.0               fastmap_1.1.0                
[137] httr_1.4.2                    survival_3.2-11              
[139] interactiveDisplayBase_1.30.0 glue_1.4.2                   
[141] png_0.1-7                     iterators_1.0.13             
[143] BiocVersion_3.13.1            bit_4.0.4                    
[145] stringi_1.6.2                 blob_1.2.1                   
[147] BiocSingular_1.8.0            AnnotationHub_3.0.0          
[149] memoise_2.0.0                 IRkernel_1.2                 
[151] irlba_2.3.3                   future.apply_1.7.0